AIM-AHEAD Project

Exploratory data analysis

Data filtering

Patient Enrollment Filtering: Claims data were filtered to include only patients with enrollment periods exceeding 12 months. This criterion was set to ensure that the data represented a sufficiently long period to analyze treatment patterns and outcomes accurately.

Linking rx_information: The rx_information data was linked to the claims data using National Drug Codes (NDCs). This linkage was essential for integrating prescription information with the claims data.

Identifying Buprenorphine OUD Cohort: The buprenorphine opioid use disorder (OUD) cohort was identified using the Universal Service Code (usc_cd) of 78430, which corresponds to ‘Drug Dependence’.

Chunk Processing Function: A custom function was written to process each chunk of data separately. This approach ensured that the dataset was handled efficiently and effectively, especially given the potentially large size of the claims data. The processed chunks were then combined into a single dataset.

Duplicate Removal: Duplicate entries were meticulously removed to prevent redundancy and ensure the accuracy of the dataset.

Joining with Demographics: The resulting data frame was joined with demographic information to create the final dataset, named ‘claims_demo’. This final dataset integrated all necessary data elements, including claims, prescriptions, and patient demographics, providing a comprehensive foundation for subsequent analyses.

The final cohort now has 3358850 claims total and 187790 individual patients.

Data linkage

Linking Demographic Data with Claims: The demographic data was linked to the claims data via patient IDs. This step ensured that each claim could be associated with the corresponding patient’s demographic information.

Linking RX Information with Claims: The rx_information data was linked to the claims data using the NDCs. This linkage was essential for associating specific medications with patient claims.

‘claims_demo’ is the final dataset.

library(dplyr)
library(data.table)
library(ggplot2) 
library(plotly)
library(usmap)
library(gridExtra)

Attaching package: ‘gridExtra’

The following object is masked from ‘package:randomForest’:

    combine

The following object is masked from ‘package:dplyr’:

    combine

filter_cohort<-function(df) {
class(df$ndc) <- "character"
class(rx_information$ndc) <- "character"
claims_rx <- df %>% 
  left_join(.,rx_information, by = "ndc")
keywords <- c("BUNAVAIL", "BUPRENORPHINE/NALOX", "BUPRENORPHINE HCL", "BUPRENORPHINE HCL/NALOXON", "BUPRENORPHINE HYDROCHLORI", "PROBUPHINE IMPLANT KIT", "SUBLOCADE", "SUBOXONE", "SUBUTEX", "ZUBSOLV") 
pattern<-paste(keywords,collapse="|")
claims_rxb <-claims_rx %>% 
  filter(grepl(pattern, product_name, ignore.case=TRUE) | grepl(pattern, generic_name, ignore.case=TRUE))
claims_rxb_12months <- claims_rxb %>% 
  filter(.$pat_id %in% patients_full_enroll)
duplicates<-apply(claims_rxb_12months,1,function(row) length(unique(row)==1))
claims_rxb_12months<-claims_rxb_12months[!duplicates, ]
claims_rxb_12mo<-claims_rxb_12months[claims_rxb_12months$usc_cd == "78340", ]
fwrite(claims_rxb_12mo, "df2_filtered")
claims_rxb_12mo2<-claims_rxb_12mo %>% select(pat_id, dayssup, quan, from_dt, to_dt, ndc, proc_cde, diagprc_ind, diag1, diag2, diag3, diag4, diag5, diag6, diag7, diag8) %>% arrange(pat_id)
fwrite(claims_rxb_12mo2, "df2_filtered_selected") 
return(head(claims_rxb_12mo2)) 
} 
filter_cohort(d2)

missingness check

missingness <- colSums(is.na(claims_demo[,c(1:20)])) / nrow(claims_demo[,c(1:20)]) * 100

missingness_df <- data.frame(
  Column = names(missingness),
  MissingPercentage = missingness
) %>%
  arrange(desc(MissingPercentage))

p <- plot_ly(
  data = missingness_df,
  x = ~MissingPercentage,
  y = ~reorder(Column, MissingPercentage),
  type = 'bar',
  orientation = 'h',
  text = ~paste("Missing: ", round(MissingPercentage, 2), "%"),
  hoverinfo = 'text',
  marker = list(color = 'steelblue')
) %>%
  layout(
    title = 'Percentage of Missing Values by Column',
    xaxis = list(title = 'Missing Percentage (%)'),
    yaxis = list(title = 'Column', automargin = TRUE, tickfont = list(size = 10)),
    margin = list(l = 200, r = 50, t = 50, b = 50)  # Adjust margins to fit labels
  )

p

missingness <- colSums(is.na(claims_demo[,-c(1:20)])) / nrow(claims_demo[,-c(1:20)]) * 100

missingness_df <- data.frame(
  Column = names(missingness),
  MissingPercentage = missingness
) %>%
  arrange(desc(MissingPercentage))

p <- plot_ly(
  data = missingness_df,
  x = ~MissingPercentage,
  y = ~reorder(Column, MissingPercentage),
  type = 'bar',
  orientation = 'h',
  text = ~paste("Missing: ", round(MissingPercentage, 2), "%"),
  hoverinfo = 'text',
  marker = list(color = 'steelblue')
) %>%
  layout(
    title = 'Percentage of Missing Values by Column',
    xaxis = list(title = 'Missing Percentage (%)'),
    yaxis = list(title = 'Column', automargin = TRUE, tickfont = list(size = 10)),
    margin = list(l = 200, r = 50, t = 50, b = 50)  # Adjust margins to fit labels
  )

p

Since diag9-diag12 are missing 100%, we could remove those columns

claims_demo<- claims_demo %>%
  select(-diag9, -diag10, -diag11, -diag12)

Since diag9-diag12 are missing 100%, we could remove those columns

Group by patient ID and dates claimed

claims_demo <- claims_demo %>% arrange(pat_id, from_dt) %>% group_by(pat_id) 

claims_demo <- claims_demo %>%
  mutate(year = year(ymd(from_dt)))

The final cohort ready for identifying episodes

head(claims_demo,100)

Summaries of Days prescribed and quantity prescribed

table(claims_demo$dayssup)

    -90     -60     -35     -33     -32     -30     -29     -28 
      2       4       1       1       1     247       2      24 
    -27     -26     -25     -24     -23     -22     -21     -20 
      1       3       1       1       1       2       6      10 
    -18     -16     -15     -14     -12     -11     -10      -9 
      2       5      26      30       1       1      21       6 
     -8      -7      -6      -5      -4      -3      -2      -1 
     16      26      19     129      18      17      33      28 
      0       1       2       3       4       5       6       7 
   7569   35993   43963   45592   40734   59788   40755  392620 
      8       9      10      11      12      13      14      15 
  63300   15387   61231   12231   19830   15489  303467  140283 
     16      17      18      19      20      21      22      23 
  17614    7371    9899    4454   35576   56336   18637   14852 
     24      25      26      27      28      29      30      31 
  12068   26315   14936   14248  391483   16661 1391960    2739 
     32      33      34      35      36      37      38      39 
   2205    1121    1818    1670     350     238     141      34 
     40      41      42      43      44      45      46      47 
    714      17     331      29      34     721      26      12 
     48      49      50      51      52      53      54      55 
     45      18     170      11      14      17      14      14 
     56      57      58      59      60      61      62      63 
    245       6       9       4    1703       9       3      12 
     64      65      66      67      68      69      70      71 
     20       8       3       4       6       1      30       2 
     72      73      75      76      77      78      80      82 
      7       3      54      13       3       3      49       1 
     83      84      85      86      87      88      89      90 
      5      73      13       1       1      13       3    3784 
     92      93      95      96      99     100     105     112 
      1       5       1       2       1      11       1       3 
    118     120     128     140     141     150     180     240 
      1       7       1       1       1       5      18       2 
    257     300     301     320     328     500     999 
      1       7       2       1       1       2       3 

dayssup_freq <- table(claims_demo$dayssup)


df_dayssup <- as.data.frame(dayssup_freq)
colnames(df_dayssup) <- c("dayssup", "frequency")

n <- nrow(df_dayssup)
part_size <- ceiling(n / 4)

df_dayssup_part1 <- df_dayssup[1:part_size, ]
df_dayssup_part2 <- df_dayssup[(part_size + 1):(2 * part_size), ]
df_dayssup_part3 <- df_dayssup[(2 * part_size + 1):(3 * part_size), ]
df_dayssup_part4 <- df_dayssup[(3 * part_size + 1):n, ]

plot1 <- ggplot(df_dayssup_part1, aes(x = as.factor(dayssup), y = frequency)) +
  geom_bar(stat = "identity", fill = "skyblue") +
  geom_text(aes(label = frequency), vjust = -0.3, size = 3.5) +
  labs(x = "Days Supply", y = "Frequency", title = "Histogram of Days Supply (Part 1)") +
  theme_minimal() +
  theme(axis.text.x = element_text(angle = 45, hjust = 1))

plot2 <- ggplot(df_dayssup_part2, aes(x = as.factor(dayssup), y = frequency)) +
  geom_bar(stat = "identity", fill = "skyblue") +
  geom_text(aes(label = frequency), vjust = -0.3, size = 3.5) +
  labs(x = "Days Supply", y = "Frequency", title = "Histogram of Days Supply (Part 2)") +
  theme_minimal() +
  theme(axis.text.x = element_text(angle = 45, hjust = 1))

plot3 <- ggplot(df_dayssup_part3, aes(x = as.factor(dayssup), y = frequency)) +
  geom_bar(stat = "identity", fill = "skyblue") +
  geom_text(aes(label = frequency), vjust = -0.3, size = 3.5) +
  labs(x = "Days Supply", y = "Frequency", title = "Histogram of Days Supply (Part 3)") +
  theme_minimal() +
  theme(axis.text.x = element_text(angle = 45, hjust = 1))

plot4 <- ggplot(df_dayssup_part4, aes(x = as.factor(dayssup), y = frequency)) +
  geom_bar(stat = "identity", fill = "skyblue") +
  geom_text(aes(label = frequency), vjust = -0.3, size = 3.5) +
  labs(x = "Days Supply", y = "Frequency", title = "Histogram of Days Supply (Part 4)") +
  theme_minimal() +
  theme(axis.text.x = element_text(angle = 45, hjust = 1))

print(plot1)

print(plot2)

print(plot3)

print(plot4)

table(claims_demo$quan)

   -90    -60    -58    -56    -28    -26    -12     -7      0 
     2      1      1      1      1      1      1      1  54831 
   0.5    0.7    0.9      1    1.1  1.111    1.3    1.5      2 
  2442      1      2  11794      3      3      2   2707  25628 
   2.5      3    3.5      4    4.5      5    5.5    5.6      6 
     1  22863      3  31202      1  26540      3      1  32881 
   6.5   6.99      7   7.03    7.5      8      9    9.5    9.9 
     1      1  86479      1      4  31446  15973      1      1 
    10   10.5     11     12   12.5     13   13.5     14     15 
 56783      9  45648  30625      4  11793      1 254479  87155 
    16  16.04     17   17.5     18     19     20     21     22 
 23352      1   7501      2  34466   4403  38616 106244  14743 
  22.5     23   23.8     24   24.5     25     26     27   27.5 
     5  16811      1  16383      1  17791   8446   6888      1 
    28     29   29.7     30     31     32     33     34     35 
189114   3517      1 345997   2176  12683   3809   4324  29121 
    36     37   37.5     38     39     40     41     42     43 
  7289   4231      1  12975   3054  15890   1593 102241   1562 
    44     45     46     47     48     49     50     51     52 
  4142 158437   3579    979   4478   9130  13049   1469   4969 
  52.5     53     54     55     56     57     58     59     60 
     4   7944   5445   3614 155650   1660   6326    649 509950 
    61     62     63     64     65     66     67     68     69 
   312   1463   9027   1249   2593   1517    753   3054    884 
    70     71     72     73     74     75     76     77     78 
 26972    241   1439    760    659  63298    543   1556   1256 
    79     80     81     82     83     84     85     86     87 
   209   3551   1063    609   1386  47881   1455    621   1660 
    88     89     90     91     92     93     94     95     96 
   737    114 271522    217    165    477     86    278    292 
    97     98     99    100    101    102    103    104    105 
    97   1543     78   3707     22   1566     16    131   6096 
   106    107    108    109    110    111    112    113    114 
    77     31    239     19    406     42   5471     76     92 
   115    116    117    118    119    120    121    122    123 
   182    260     23     48     31  56357     10     28     31 
   124    125    126    127    128    129    130    131    132 
   120     65    254     17     60     11    231      2     38 
   133    134    135    136    137    138    139    140    142 
     2     34    967     56      3     23      2    489     10 
   143    144    145    146    147    148    149    150    151 
     4     71     42      8     29     12      1   6053      1 
   152    153    154    155    156    157    158    159    160 
     7      5     37     10     11      3      4      5     64 
   162    163    164    165    166    167    168    170    172 
    12      1      9    113      7      1    399     31     10 
   173    174    175    176    177    178    179    180    181 
     1     14      7      5      4      3      1   7300      1 
   182    184    185    186    187    189    190    192    195 
    11      3      5      7      1      5     15     11     36 
   196    198    200    201    203    204    205    206    208 
   267      1    108      1      7     11      2      1      3 
   209    210    211    212    215    216    217    218    220 
     6    726      1      1      2      5      1      1     38 
   221    224    225    228    230    231    232    240    243 
     1    134    129      1      9      1     13   1757      2 
   244    245    246    250    252    255    256    257    260 
     2      1      1      8     17      3      1      1      8 
   261    264    265    268    269    270    275    280    290 
     1      1      1      1      1   1267      4     37     11 
   294    300    304    306    315    320    322    324    325 
     2    321      1      1      7      3      1      2      2 
   330    333    336    340    341    348    350    354    360 
    25      1     15      1      1      3      4      4    577 
   370    375    380    384    390    400    405    410    420 
     3      3     22      1     20      4      6      1     29 
   440    442    448    450    464    476    480    500    504 
     1      1      1     82      2      1     60      1      1 
   510    520    540    543    550    560    570    580    600 
     3      3     44      2      1     23      1      3    162 
   620    630    660    680    690    700    720    750    780 
     1      2      1      1      1      3      6     13      3 
   800    810    840    870    890    900    930    960    980 
     2     14     12      1      1     62      1      1      5 
  1000   1050   1060   1080   1100   1120   1200   1230   1800 
     3      5      2      1      1      5     20      1      2 
  2000   3000   5000   5829   7000   8000   9000  11000  12000 
    10     10      8      1     12      4      2      6     10 
 13000  14000  16000  21000  22000  23000  24000  25000  26000 
     2      5      6     20      9     12      1     18      1 
 27000  32000  33000  35000  38000  40000  41000  42000  43000 
     2     10      4      2      4      6      2     13      1 
 44000  45000  46000  50000  53000  56000  62000  63000  66000 
     2     63      3      2      6      1      1      3      2 
270000 
     3 

quantities_freq <- table(claims_demo$quan) 
df_quantities <- as.data.frame(quantities_freq)
colnames(df_quantities) <- c("quantities", "frequency")

n <- nrow(df_quantities)
part_size <- ceiling(n / 8)

df_quantities_part1 <- df_quantities[1:part_size, ]
df_quantities_part2 <- df_quantities[(part_size + 1):(2 * part_size), ]
df_quantities_part3 <- df_quantities[(2 * part_size + 1):(3 * part_size), ]
df_quantities_part4 <- df_quantities[(3 * part_size + 1):(4 * part_size), ]
df_quantities_part5 <- df_quantities[(4 * part_size + 1):(5 * part_size), ]
df_quantities_part6 <- df_quantities[(5 * part_size + 1):(6 * part_size), ]
df_quantities_part7 <- df_quantities[(6 * part_size + 1):(7 * part_size), ]
df_quantities_part8 <- df_quantities[(7 * part_size + 1):n, ]

plot1 <- ggplot(df_quantities_part1, aes(x = as.factor(quantities), y = frequency)) +
  geom_bar(stat = "identity", fill = "skyblue") +
  geom_text(aes(label = frequency), vjust = -0.3, size = 3.5) +
  labs(x = "Quantities", y = "Frequency", title = "Histogram of Quantities (Part 1)") +
  theme_minimal() +
  theme(axis.text.x = element_text(angle = 45, hjust = 1))

plot2 <- ggplot(df_quantities_part2, aes(x = as.factor(quantities), y = frequency)) +
  geom_bar(stat = "identity", fill = "skyblue") +
  geom_text(aes(label = frequency), vjust = -0.3, size = 3.5) +
  labs(x = "Quantities", y = "Frequency", title = "Histogram of Quantities (Part 2)") +
  theme_minimal() +
  theme(axis.text.x = element_text(angle = 45, hjust = 1))

plot3 <- ggplot(df_quantities_part3, aes(x = as.factor(quantities), y = frequency)) +
  geom_bar(stat = "identity", fill = "skyblue") +
  geom_text(aes(label = frequency), vjust = -0.3, size = 3.5) +
  labs(x = "Quantities", y = "Frequency", title = "Histogram of Quantities (Part 3)") +
  theme_minimal() +
  theme(axis.text.x = element_text(angle = 45, hjust = 1))

plot4 <- ggplot(df_quantities_part4, aes(x = as.factor(quantities), y = frequency)) +
  geom_bar(stat = "identity", fill = "skyblue") +
  geom_text(aes(label = frequency), vjust = -0.3, size = 3.5) +
  labs(x = "Quantities", y = "Frequency", title = "Histogram of Quantities (Part 4)") +
  theme_minimal() +
  theme(axis.text.x = element_text(angle = 45, hjust = 1))

plot5 <- ggplot(df_quantities_part5, aes(x = as.factor(quantities), y = frequency)) +
  geom_bar(stat = "identity", fill = "skyblue") +
  geom_text(aes(label = frequency), vjust = -0.3, size = 3.5) +
  labs(x = "Quantities", y = "Frequency", title = "Histogram of Quantities (Part 5)") +
  theme_minimal() +
  theme(axis.text.x = element_text(angle = 45, hjust = 1))

plot6 <- ggplot(df_quantities_part6, aes(x = as.factor(quantities), y = frequency)) +
  geom_bar(stat = "identity", fill = "skyblue") +
  geom_text(aes(label = frequency), vjust = -0.3, size = 3.5) +
  labs(x = "Quantities", y = "Frequency", title = "Histogram of Quantities (Part 6)") +
  theme_minimal() +
  theme(axis.text.x = element_text(angle = 45, hjust = 1))

plot7 <- ggplot(df_quantities_part7, aes(x = as.factor(quantities), y = frequency)) +
  geom_bar(stat = "identity", fill = "skyblue") +
  geom_text(aes(label = frequency), vjust = -0.3, size = 3.5) +
  labs(x = "Quantities", y = "Frequency", title = "Histogram of Quantities (Part 7)") +
  theme_minimal() +
  theme(axis.text.x = element_text(angle = 45, hjust = 1))

plot8 <- ggplot(df_quantities_part8, aes(x = as.factor(quantities), y = frequency)) +
  geom_bar(stat = "identity", fill = "skyblue") +
  geom_text(aes(label = frequency), vjust = -0.3, size = 3.5) +
  labs(x = "Quantities", y = "Frequency", title = "Histogram of Quantities (Part 8)") +
  theme_minimal() +
  theme(axis.text.x = element_text(angle = 45, hjust = 1))

print(plot1)

print(plot2)

print(plot3)

print(plot4)

print(plot5)

print(plot6)

print(plot7)

print(plot8)

Summarize buprenorphine prescriptions by year

prescriptions_by_year <- claims_demo %>%
  group_by(year) %>%
  summarise(count = n())

#1. Bar Graph of Counts by Year
# Plot the bar graph
ggplot(prescriptions_by_year, aes(x = year, y = count)) +
  geom_bar(stat = "identity") +
  labs(title = "Number of Prescriptions by Year",
       x = "Year",
       y = "Number of Prescriptions") +
  theme_minimal() +
  theme(axis.text.x = element_text(angle = 45, hjust = 1))

From 2006 to 2013, there was a steady increase in the number of prescriptions, peaking in 2013 and 2014 with over 300,000 prescriptions each year.

This rise suggests a growing trend in prescription rates or improved access to healthcare during these years.

However, following 2014, there is a notable decline in the number of prescriptions, reaching its lowest point around 2019. This decline may be attributed to changes in healthcare policies.

In the subsequent years, there is a slight recovery in 2021 and 2022, indicating a resurgence in prescription rates, but this is followed by a sharp decrease in 2023. This could be influenced by external factors such as changes in healthcare regulations, economic conditions, or the impact of global events like the COVID-19 pandemic.

Summarize buprenorphine prescriptions by year, Stratified by Sex

Filter the data to include only entries with sex “F” and “M”

fig1_filtered <- claims_demo %>%
  filter(der_sex %in% c("F", "M"))

# Summarize the data by year and sex
prescriptions_by_year_sex <- fig1_filtered %>%
  group_by(year, der_sex) %>%
  summarise(count = n())

`summarise()` has grouped output by 'year'. You can override using the `.groups` argument.

# Plot the bar graph
ggplot(prescriptions_by_year_sex, aes(x = year, y = count, fill = der_sex)) +
  geom_bar(stat = "identity", position = "dodge") +
  labs(title = "Number of Prescriptions by Year, Stratified by Sex",
       x = "Year",
       y = "Number of Prescriptions",
       fill = "Sex") +
  theme_minimal() +
  theme(axis.text.x = element_text(angle = 45, hjust = 1))

The second bar graph shows the number of prescriptions by year, stratified by sex (F and M). It highlights several important trends and differences between male and female prescription rates. Overall, the number of prescriptions for both sexes follows a similar pattern, with an increase from 2006, peaking around 2013 and 2014, followed by a decline.

Throughout the years, males consistently have higher prescription counts compared to females. The peak years of 2013 and 2014 show the most significant differences, with males having substantially more prescriptions than females. After 2014, the number of prescriptions for both sexes declines, reaching a low point around 2019.

This graph indicates that while the overall prescription trends are similar for both males and females, males consistently receive more prescriptions than females across all years.

create a map of prescriptions by US state

fig2<-claims_demo %>% select(pat_state, year) 
write.csv(fig2, "fig2.csv")  

prescriptions_by_state <- fig2 %>%
  group_by(pat_state) %>%
  summarise(total_prescriptions = n())

head(prescriptions_by_state)
prescriptions_by_state <- as.data.frame(prescriptions_by_state)
prescriptions_by_state<-prescriptions_by_state[-1,]
colnames(prescriptions_by_state)[1]<-"state"

plot_usmap(data = prescriptions_by_state, values = "total_prescriptions", color = "white") +
  scale_fill_continuous(name = "Total Prescriptions", label = scales::comma) +
  labs(title = "Total Number of Prescriptions by US State") +
  theme(legend.position = "right")

create a map of Total Prescriptions by US State

state_data <- merge(prescriptions_by_state, state_centers, by.x = "state", by.y = "state")

# Total prescriptions plot using plotly
total_prescriptions_plot <- plot_usmap(data = prescriptions_by_state, values = "total_prescriptions", color = "white") +
  geom_text(data = state_data, aes(x = x, y = y, label = state), size = 3, hjust = 0.5, color = "black") +
  scale_fill_gradientn(name = "Total Prescriptions", colors = c("blue", "green", "yellow", "red"), labels = scales::comma, 
                       breaks = seq(0, max(prescriptions_by_state$total_prescriptions), length.out = 6)) +
  labs(title = "Total Number of Prescriptions by US State") +
  theme(legend.position = "right")

# Convert to plotly object
ggplotly(total_prescriptions_plot)

Create a Shiny App

ui <- fluidPage(
  titlePanel("Number of Prescriptions by US State"),
  tabsetPanel(
    tabPanel("Total",
             plotOutput("staticMap")
    ),
    tabPanel("Yearly",
             sidebarLayout(
               sidebarPanel(
                 selectInput("year", "Select Year:", choices = unique(prescriptions_by_state_year$year))
               ),
               mainPanel(
                 plotlyOutput("usMap")
               )
             )
    ),
    tabPanel("Product Trends",
             plotOutput("productTrendGraph")
    )
  )
)

server <- function(input, output) {
  output$staticMap <- renderPlot({
    state_data <- merge(prescriptions_by_state, state_centers, by.x = "state", by.y = "state")
    
    plot_usmap(data = prescriptions_by_state, values = "total_prescriptions", color = "white") +
      geom_text(data = state_data, aes(x = x, y = y, label = state), size = 3, hjust = 0.5, color = "black") +
      scale_fill_gradientn(name = "Total Prescriptions", colors = c("blue", "green", "yellow", "red"), labels = scales::comma, 
                           breaks = seq(0, max(prescriptions_by_state$total_prescriptions), length.out = 6)) +
      labs(title = "Total Number of Prescriptions by US State") +
      theme(legend.position = "right")
  })
  
  output$usMap <- renderPlotly({
    data_year <- prescriptions_by_state_year %>%
      filter(year == input$year)
    
    state_data_year <- merge(data_year, state_centers, by = "state")
    
    p <- plot_usmap(data = data_year, values = "total_prescriptions", color = "white") +
      geom_text(data = state_data_year, aes(x = x, y = y, label = state), size = 3, hjust = 0.5, color = "black") +
      scale_fill_gradientn(name = paste("Prescriptions in", input$year), colors = c("blue", "green", "yellow", "red"), labels = scales::comma, 
                           breaks = seq(0, max(data_year$total_prescriptions), length.out = 6)) +
      labs(title = paste("Number of Buprenorphine OUD Prescriptions by US State in", input$year)) +
      theme(legend.position = "right")
    
    ggplotly(p)
  })
  
  output$productTrendGraph <- renderPlot({
    ggplot(product_data, aes(x = year, y = percentage, color = product_name)) +
      geom_line(size = 1.2) +
      scale_color_viridis(discrete = TRUE) + 
      labs(title = "Percentage of Each Product Over Time (2006-2023)",
           x = "Year",
           y = "Percentage of Total Prescriptions",
           color = "Product") +
      theme_minimal() +
      theme(legend.position = "right", legend.title = element_text(size = 10),
            legend.text = element_text(size = 8))
  })
}

shinyApp(ui = ui, server = server)

Create a temporal graph of Product trends Over Time

product_trend_plot <- ggplot(product_data, aes(x = year, y = percentage, color = product_name)) +
  geom_line(size = 1.2) +
  scale_color_brewer(palette = "Set3") + 
  labs(title = "Percentage of Each Product Over Time (2006-2023)",
       x = "Year",
       y = "Percentage of Total Prescriptions",
       color = "Product") +
  theme_minimal() +
  theme(legend.position = "right", legend.title = element_text(size = 10),
        legend.text = element_text(size = 8))

# Convert to plotly object
ggplotly(product_trend_plot)

240923*0.572 ## discussion of product Trends 240923-137808 100-57.2 866904/240923

SUBOXONE Trend: There is a significant increase in the percentage of SUBOXONE prescriptions starting around 2010, reaching a peak around 2015-2016, followed by a sharp decline.

Literature Support: SUBOXONE, a combination of buprenorphine and naloxone, was approved by the FDA in 2002 for the treatment of opioid dependence. The increase in its prescription aligns with its effectiveness in treating opioid use disorder (OUD) and its expanded use in clinical practice. The decline in recent years could be attributed to the introduction of new buprenorphine products and the implementation of generic versions, which tend to be less expensive and more accessible.

BUPRENORPHINE HCL Trend: Shows a significant increase starting around 2018. Literature Support: Buprenorphine HCL has been used for both pain management and OUD. The increase in recent years can be attributed to the broader acceptance and understanding of buprenorphine’s efficacy in treating OUD. It may also reflect the FDA’s efforts to expand access to buprenorphine as part of the response to the opioid crisis.

BUPRENORPHINE HCL/NALOXONE Trend: Steady increase in percentage over time after 2013. Literature Support: This trend supports the literature that emphasizes the effectiveness of combination therapies (buprenorphine/naloxone) in reducing the potential for misuse compared to buprenorphine alone. The increase in prescriptions also coincides with the FDA’s encouragement of such combinations for safety reasons.

SUBUTEX Trend: High percentage initially, followed by a decline around 2010 and remaining low thereafter. Literature Support: SUBUTEX, which contains only buprenorphine, was commonly used before the introduction of combination products. Its decline is supported by FDA recommendations favoring combination products like SUBOXONE to reduce misuse potential.

SUBLOCADE Trend: Emerges around 2018 with a rising trend. Literature Support: SUBLOCADE, an extended-release injectable form of buprenorphine, was approved by the FDA in 2017. Its increasing prescription rate reflects the growing acceptance and utilization of long-acting formulations for better compliance and sustained recovery in OUD treatment.

PROBUPHINE IMPLANT KIT Trend: Shows a very low and relatively flat trend. Literature Support: PROBUPHINE, approved in 2016, is an implantable form of buprenorphine designed to provide continuous delivery over six months. Its limited use might be due to the invasive nature of implantation and the preference for less invasive administration routes.

ZUBSOLV Trend: Shows a minor presence with some increase over time. Literature Support: ZUBSOLV is another buprenorphine/naloxone combination product with a different formulation and bioavailability. Its introduction provided an alternative to SUBOXONE, and I believe is highly reliant on prescriber habits.

BUNAVAIL Trend: Minor presence with slight fluctuation. Literature Support: BUNAVAIL, a buccal film formulation of buprenorphine/naloxone, was introduced as an alternative to sublingual formulations. Its relatively minor uptake could be due to competition with other established combination products.

Trends Not Fully Supported by Literature

Sharp Decline in SUBOXONE: The rapid decline in SUBOXONE prescriptions post-2015 is sharper than might be expected solely due to the introduction of generics and new formulations. This could indicate other factors at play, such as changes in insurance coverage, prescriber preferences, or regulatory actions targeting specific products.

Relatively Flat Trend for PROBUPHINE: Despite its innovative implant approach, PROBUPHINE’s limited adoption might not fully align with its potential benefits as highlighted in clinical studies. Factors such as cost, patient and provider acceptance, and procedural complexity likely play a role.

summarizing the final cohort dataset

rectype_summary <- summary(claims_demo$rectype)
dosage_form_nm_summary <- summary(claims_demo$dosage_form_nm)
route_summary <- summary(claims_demo$route)
strength_summary <- summary(claims_demo$strength)

rectype_plot <- ggplot(claims_demo, aes(x = factor(rectype))) +
  geom_bar(fill = 'skyblue', color = 'black') +
  labs(title = 'Distribution of Rectype', x = 'Rectype', y = 'Frequency') +
  theme_minimal()

dosage_form_nm_plot <- ggplot(claims_demo, aes(x = factor(dosage_form_nm))) +
  geom_bar(fill = 'salmon', color = 'black') +
  labs(title = 'Distribution of Dosage Form Name', x = 'Dosage Form Name', y = 'Frequency') +
  theme_minimal() +
  theme(axis.text.x = element_text(angle = 45, hjust = 1))


route_plot <- ggplot(claims_demo, aes(x = factor(route))) +
  geom_bar(fill = 'lightgreen', color = 'black') +
  labs(title = 'Distribution of Route', x = 'Route', y = 'Frequency') +
  theme_minimal() +
  theme(axis.text.x = element_text(angle = 45, hjust = 1))


strength_plot <- ggplot(claims_demo, aes(x = factor(strength))) +
  geom_bar(fill = 'lightblue', color = 'black') +
  labs(title = 'Distribution of Strength', x = 'Strength', y = 'Frequency') +
  theme_minimal() +
  theme(axis.text.x = element_text(angle = 45, hjust = 1))

rectype_plot 

dosage_form_nm_plot 

route_plot 

strength_plot 

colnames(claims_demo)

selected_data <- claims_demo %>% select(pat_id, der_sex, pat_region, pat_state, yob_group, zip3_group)

Error in `select()`:
! Can't subset columns that don't exist.
✖ Column `yob_group` doesn't exist.
Backtrace:
 1. claims_demo %>% ...
 3. dplyr:::select.data.frame(., pat_id, der_sex, pat_region, pat_state, yob_group, zip3_group)